# Replication code for Claudia J. Kim and Taylor C. Boas, "Activist Disconnect: Social Movements, Public Opinion, and U.S. Military Bases in East Asia," Armed Forces and Society.

# Analysis conducted in R 3.6.0 on MacOS 10.13.6

# NOTE: This file reproduces Appendix Table 1 and Appendix Figures 3-6. 

# Set working directory as appropriate
# setwd('~/Dropbox/us_bases_experiment/replication/')

# Clean desktop and load packages. Please make sure all necessary packages are installed.

rm(list=ls(all=T))

library(readxl)
library(Hmisc)
library(foreign)

# Load cleaned dataset, census data, and World Values Survey data for Japan and South Korea

load('basedata_cleaned.RData')

daegu_age<-as.data.frame(read_excel('age_census_new.xlsx',sheet=1,skip=1))
gyeonggi_age<-as.data.frame(read_excel('age_census_new.xlsx',sheet=2,skip=1))
okinawa_age<-as.data.frame(read_excel('age_census_new.xlsx',sheet=3,skip=1))
kanagawa_age<-as.data.frame(read_excel('age_census_new.xlsx',sheet=4,skip=1))

korea<-read.spss('WV6_Data_South_Korea_2010_Spss_v20180912.sav')
japan<-read.spss('WV6_Data_Japan_2010_Spss_v20180912.sav')

##################
# Appendix Table 1
##################

# Population age and sex distribution

daegu_age<-daegu_age[!is.na(daegu_age$age),1:5]
names(daegu_age)<-c('age','men','men_pct','women','women_pct')
daegu_age$total<-daegu_age$men + daegu_age$women
daegu_age$total_pct<-100* daegu_age$total / sum(daegu_age$total)
daegu_age$cumpct<-cumsum(daegu_age$total_pct)

gyeonggi_age<-gyeonggi_age[!is.na(gyeonggi_age$age),1:5]
names(gyeonggi_age)<-c('age','men','men_pct','women','women_pct')
gyeonggi_age$total<-gyeonggi_age$men + gyeonggi_age$women
gyeonggi_age$total_pct<-100* gyeonggi_age$total / sum(gyeonggi_age$total)
gyeonggi_age$cumpct<-cumsum(gyeonggi_age$total_pct)

okinawa_age<-okinawa_age[!is.na(okinawa_age$age),1:5]
names(okinawa_age)<-c('age','men','men_pct','women','women_pct')
okinawa_age$total<-okinawa_age$men + okinawa_age$women
okinawa_age$total_pct<-100* okinawa_age$total / sum(okinawa_age$total)
okinawa_age$cumpct<-cumsum(okinawa_age$total_pct)

kanagawa_age<-kanagawa_age[!is.na(kanagawa_age$age),1:5]
names(kanagawa_age)<-c('age','men','men_pct','women','women_pct')
kanagawa_age$total<-kanagawa_age$men + kanagawa_age$women
kanagawa_age$total_pct<-100* kanagawa_age$total / sum(kanagawa_age$total)
kanagawa_age$cumpct<-cumsum(kanagawa_age$total_pct)

# Median ages by region

daegu_med_age<-as.numeric(daegu_age$age[which(daegu_age$cumpct >= 50)[1]])
gyeonggi_med_age<-as.numeric(gyeonggi_age$age[which(gyeonggi_age$cumpct >= 50)[1]])
okinawa_med_age<-as.numeric(okinawa_age$age[which(okinawa_age$cumpct >= 50)[1]])
kanagawa_med_age<-as.numeric(kanagawa_age$age[which(kanagawa_age$cumpct >= 50)[1]])

# Combo median age, weighting each region equally

combo_age<-data.frame(age=daegu_age$age,cumpct=apply(cbind(daegu_age$cumpct, gyeonggi_age$cumpct, okinawa_age$cumpct, kanagawa_age$cumpct),1,mean),stringsAsFactors=F)
combo_med_age<-as.numeric(combo_age$age[which(combo_age$cumpct >= 50)[1]])

# Percent male by region

daegu_pct_male<-100*sum(daegu_age$men) / sum(c(daegu_age$men, daegu_age$women))
gyeonggi_pct_male<-100*sum(gyeonggi_age$men) / sum(c(gyeonggi_age$men, gyeonggi_age$women))
okinawa_pct_male<-100*sum(okinawa_age$men) / sum(c(okinawa_age$men, okinawa_age$women))
kanagawa_pct_male<-100*sum(kanagawa_age$men) / sum(c(kanagawa_age$men, kanagawa_age$women))

# Combo percent male, weighting each region equally

combo_pct_male<-mean(daegu_pct_male, gyeonggi_pct_male, okinawa_pct_male, kanagawa_pct_male)

# Population by region and combined

daegu_pop<-sum(c(daegu_age$men, daegu_age$women))
gyeonggi_pop<-sum(c(gyeonggi_age$men, gyeonggi_age$women))
okinawa_pop<-sum(c(okinawa_age$men, okinawa_age$women))
kanagawa_pop<-sum(c(kanagawa_age$men, kanagawa_age$women))
combo_pop<-daegu_pop+ gyeonggi_pop+ okinawa_pop+ kanagawa_pop

# Sample age and sex distribution

basedata$Age<-basedata$Q4
basedata$Male<-basedata$Q33==1
basedata$Region<-basedata$Q5

reg_weight<-data.frame((1/4)/prop.table(table(basedata$Region)))
names(reg_weight)<-c('Region','reg_weight')
basedata<-merge(basedata,reg_weight)

daegu_med_age_sample<-median(basedata$Age[basedata$Region=='Daegu'])
gyeonggi_med_age_sample<-median(basedata$Age[basedata$Region=='Gyeonggi'])
kanagawa_med_age_sample<-median(basedata$Age[basedata$Region=='Kanagawa'])
okinawa_med_age_sample<-median(basedata$Age[basedata$Region=='Okinawa'])

combo_med_age_sample<-wtd.quantile(basedata$Age,weights=basedata$reg_weight)[' 50%']

daegu_pct_male_sample<-100*tapply(basedata$Male,basedata$Region,mean)['Daegu']
gyeonggi_pct_male_sample<-100*tapply(basedata$Male,basedata$Region,mean)['Gyeonggi']
kanagawa_pct_male_sample<-100*tapply(basedata$Male,basedata$Region,mean)['Kanagawa']
okinawa_pct_male_sample<-100*tapply(basedata$Male,basedata$Region,mean)['Okinawa']
combo_pct_male_sample<-mean(100*tapply(basedata$Male,basedata$Region,mean))

# Appendix Table 1

rep.table<-round(matrix(c(daegu_pop, gyeonggi_pop, kanagawa_pop, okinawa_pop, combo_pop,
	table(basedata$Region), nrow(basedata),
	daegu_pct_male, gyeonggi_pct_male, kanagawa_pct_male, okinawa_pct_male, combo_pct_male,
	daegu_pct_male_sample, gyeonggi_pct_male_sample, kanagawa_pct_male_sample, okinawa_pct_male_sample, combo_pct_male_sample,
	daegu_med_age, gyeonggi_med_age, kanagawa_med_age, okinawa_med_age, combo_med_age,
	daegu_med_age_sample, gyeonggi_med_age_sample, kanagawa_med_age_sample, okinawa_med_age_sample, combo_med_age_sample),ncol=5,byrow=T))
colnames(rep.table)<-c('Daegu','Gyeonggi','Kanagawa','Okinawa','Combined')
rownames(rep.table)<-rep(c('Population','Sample'),3)
rep.table<-apply(rep.table,2,prettyNum,big.mark=',')

rep.table.note<-"Census figures are based on the 18-and-older population in each country's 2015 Census. Combined statistics weight each region equally."
	
rep.table.latex<-latex(rep.table,file='rep_table.tex', rowlabel = '', collabel.just=rep('r',ncol(rep.table)), col.just=rep('r',ncol(rep.table)), cgroup=c('South Korea','Japan',''),n.cgroup=c(2,2,1),rgroup=c('\\textit{N}','Percent Male','Median Age'), n.rgroup=c(2,2,2), caption = 'Online Sample versus Census', insert.bottom= rep.table.note, booktabs = F, ctable = T, where = "htp")

#################################
# Appendix Figures 3, 4, 5, and 6
#################################

basedata$Ideology<-basedata$Q13
basedata$petition<-basedata$Q8
basedata$boycott<-basedata$Q9
basedata$demonstration<-basedata$Q10

korea$ideology<-ifelse(unclass(korea$V95) %in% 1:5, NA, unclass(korea$V95) - 5)
korea$petition<-ifelse(unclass(korea$V85) %in% 1:5, NA, unclass(korea$V85) - 5)
korea$boycott<-ifelse(unclass(korea$V86) %in% 1:5, NA, unclass(korea$V86) - 5)
korea$demonstration<-ifelse(unclass(korea$V87) %in% 1:5, NA, unclass(korea$V87) - 5)

japan$ideology<-ifelse(unclass(japan$V95) %in% 1:5, NA, unclass(japan$V95) - 5)
japan$petition<-ifelse(unclass(japan$V85) %in% 1:5, NA, unclass(japan$V85) - 5)
japan$boycott<-ifelse(unclass(japan$V86) %in% 1:5, NA, unclass(japan$V86) - 5)
japan$demonstration<-ifelse(unclass(japan$V87) %in% 1:5, NA, unclass(japan$V87) - 5)

ideol.daegu.wvs<-mean(korea$ideology[korea$V256 == 'KR: Taegu / Daegu'])
ideol.daegu<-mean(basedata$Ideology[basedata$Region=='Daegu'])
ideol.gyeonggi.wvs<-mean(korea$ideology[korea$V256 == 'KR: Kyeonggi / Gyeonggi Do'])
ideol.gyeonggi<-mean(basedata$Ideology[basedata$Region=='Gyeonggi'])

ideol.minamikanto<-mean(japan$ideology[japan$V256 == 'JP: Minami-Kanto region'], na.rm=T)
ideol.kanagawa<-mean(basedata$Ideology[basedata$Region=='Kanagawa'])
ideol.kyushu<-mean(japan$ideology[japan$V256 == 'JP: Kyushu region'], na.rm=T)
ideol.okinawa<-mean(basedata$Ideology[basedata$Region=='Okinawa'])

petition.daegu.wvs<-mean(korea$petition[korea$V256 == 'KR: Taegu / Daegu'], na.rm=T)
petition.daegu<-mean(basedata$petition[basedata$Region=='Daegu'])
petition.gyeonggi.wvs<-mean(korea$petition[korea$V256 == 'KR: Kyeonggi / Gyeonggi Do'], na.rm=T)
petition.gyeonggi<-mean(basedata$petition[basedata$Region=='Gyeonggi'])

petition.minamikanto<-mean(japan$petition[japan$V256 == 'JP: Minami-Kanto region'], na.rm=T)
petition.kanagawa<-mean(basedata$petition[basedata$Region=='Kanagawa'])
petition.kyushu<-mean(japan$petition[japan$V256 == 'JP: Kyushu region'], na.rm=T)
petition.okinawa<-mean(basedata$petition[basedata$Region=='Okinawa'])

boycott.daegu.wvs<-mean(korea$boycott[korea$V256 == 'KR: Taegu / Daegu'], na.rm=T)
boycott.daegu<-mean(basedata$boycott[basedata$Region=='Daegu'])
boycott.gyeonggi.wvs<-mean(korea$boycott[korea$V256 == 'KR: Kyeonggi / Gyeonggi Do'], na.rm=T)
boycott.gyeonggi<-mean(basedata$boycott[basedata$Region=='Gyeonggi'])

boycott.minamikanto<-mean(japan$boycott[japan$V256 == 'JP: Minami-Kanto region'], na.rm=T)
boycott.kanagawa<-mean(basedata$boycott[basedata$Region=='Kanagawa'])
boycott.kyushu<-mean(japan$boycott[japan$V256 == 'JP: Kyushu region'], na.rm=T)
boycott.okinawa<-mean(basedata$boycott[basedata$Region=='Okinawa'])

demonstration.daegu.wvs<-mean(korea$demonstration[korea$V256 == 'KR: Taegu / Daegu'], na.rm=T)
demonstration.daegu<-mean(basedata$demonstration[basedata$Region=='Daegu'])
demonstration.gyeonggi.wvs<-mean(korea$demonstration[korea$V256 == 'KR: Kyeonggi / Gyeonggi Do'], na.rm=T)
demonstration.gyeonggi<-mean(basedata$demonstration[basedata$Region=='Gyeonggi'])

demonstration.minamikanto<-mean(japan$demonstration[japan$V256 == 'JP: Minami-Kanto region'], na.rm=T)
demonstration.kanagawa<-mean(basedata$demonstration[basedata$Region=='Kanagawa'])
demonstration.kyushu<-mean(japan$demonstration[japan$V256 == 'JP: Kyushu region'], na.rm=T)
demonstration.okinawa<-mean(basedata$demonstration[basedata$Region=='Okinawa'])

# Appendix Figure 3

pdf(file='ideology.pdf',width=9,height=6.5)
par(mfrow=c(2,4))
barplot(table(korea$ideology[korea$V256 == 'KR: Taegu / Daegu']), main=paste0('WVS: Daegu\n(mean = ', round(ideol.daegu.wvs,1), ')'), ylab='Frequency',xlab='Ideology (Left-Right)')
barplot(table(basedata$Ideology[basedata$Region=='Daegu']), main=paste0('Online Sample: Daegu\n(mean = ', round(ideol.daegu,1), ')'), ylab='Frequency',xlab='Ideology (Left-Right)')
barplot(table(korea$ideology[korea$V256 == 'KR: Kyeonggi / Gyeonggi Do']), main=paste0('WVS: Gyeonggi\n(mean = ', round(ideol.gyeonggi.wvs,1), ')'), ylab='Frequency',xlab='Ideology (Left-Right)')
barplot(table(basedata$Ideology[basedata$Region=='Gyeonggi']), main=paste0('Online Sample: Gyeonggi\n(mean = ', round(ideol.gyeonggi,1), ')'), ylab='Frequency',xlab='Ideology (Left-Right)')

barplot(table(japan$ideology[japan$V256 == 'JP: Minami-Kanto region']), main=paste0('WVS: Minami-Kanto Region\n(mean = ', round(ideol.minamikanto,1), ')'), ylab='Frequency',xlab='Ideology (Left-Right)')
barplot(table(basedata$Ideology[basedata$Region=='Kanagawa']), main=paste0('Online Sample: Kanagawa\n(mean = ', round(ideol.kanagawa,1), ')'), ylab='Frequency',xlab='Ideology (Left-Right)')
barplot(table(japan$ideology[japan$V256 == 'JP: Kyushu region']), main=paste0('WVS: Kyushu Region\n(mean = ', round(ideol.kyushu,1), ')'), ylab='Frequency',xlab='Ideology (Left-Right)')
barplot(table(basedata$Ideology[basedata$Region=='Okinawa']), main=paste0('Online Sample: Okinawa\n(mean = ', round(ideol.okinawa,1), ')'), ylab='Frequency',xlab='Ideology (Left-Right)')
dev.off()

# Appendix Figure 4

pdf(file='petition.pdf',width=9,height=6.5)
par(mfrow=c(2,4))
barplot(table(korea$petition[korea$V256 == 'KR: Taegu / Daegu']), main=paste0('WVS: Daegu\n(mean = ', round(petition.daegu.wvs,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(basedata$petition[basedata$Region=='Daegu']), main=paste0('Online Sample: Daegu\n(mean = ', round(petition.daegu,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(korea$petition[korea$V256 == 'KR: Kyeonggi / Gyeonggi Do']), main=paste0('WVS: Gyeonggi\n(mean = ', round(petition.gyeonggi.wvs,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(basedata$petition[basedata$Region=='Gyeonggi']), main=paste0('Online Sample: Gyeonggi\n(mean = ', round(petition.gyeonggi,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))

barplot(table(japan$petition[japan$V256 == 'JP: Minami-Kanto region']), main=paste0('WVS: Minami-Kanto Region\n(mean = ', round(petition.minamikanto,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(basedata$petition[basedata$Region=='Kanagawa']), main=paste0('Online Sample: Kanagawa\n(mean = ', round(petition.kanagawa,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(japan$petition[japan$V256 == 'JP: Kyushu region']), main=paste0('WVS: Kyushu Region\n(mean = ', round(petition.kyushu,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(basedata$petition[basedata$Region=='Okinawa']), main=paste0('Online Sample: Okinawa\n(mean = ', round(petition.okinawa,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
dev.off()

# Appendix Figure 5

pdf(file='boycott.pdf',width=9,height=6.5)
par(mfrow=c(2,4))
barplot(table(korea$boycott[korea$V256 == 'KR: Taegu / Daegu']), main=paste0('WVS: Daegu\n(mean = ', round(boycott.daegu.wvs,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(basedata$boycott[basedata$Region=='Daegu']), main=paste0('Online Sample: Daegu\n(mean = ', round(boycott.daegu,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(korea$boycott[korea$V256 == 'KR: Kyeonggi / Gyeonggi Do']), main=paste0('WVS: Gyeonggi\n(mean = ', round(boycott.gyeonggi.wvs,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(basedata$boycott[basedata$Region=='Gyeonggi']), main=paste0('Online Sample: Gyeonggi\n(mean = ', round(boycott.gyeonggi,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))

barplot(table(japan$boycott[japan$V256 == 'JP: Minami-Kanto region']), main=paste0('WVS: Minami-Kanto Region\n(mean = ', round(boycott.minamikanto,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(basedata$boycott[basedata$Region=='Kanagawa']), main=paste0('Online Sample: Kanagawa\n(mean = ', round(boycott.kanagawa,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(japan$boycott[japan$V256 == 'JP: Kyushu region']), main=paste0('WVS: Kyushu Region\n(mean = ', round(boycott.kyushu,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(basedata$boycott[basedata$Region=='Okinawa']), main=paste0('Online Sample: Okinawa\n(mean = ', round(boycott.okinawa,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
dev.off()

# Appendix Figure 6

pdf(file='demonstration.pdf',width=9,height=6.5)
par(mfrow=c(2,4))
barplot(table(korea$demonstration[korea$V256 == 'KR: Taegu / Daegu']), main=paste0('WVS: Daegu\n(mean = ', round(demonstration.daegu.wvs,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(basedata$demonstration[basedata$Region=='Daegu']), main=paste0('Online Sample: Daegu\n(mean = ', round(demonstration.daegu,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(korea$demonstration[korea$V256 == 'KR: Kyeonggi / Gyeonggi Do']), main=paste0('WVS: Gyeonggi\n(mean = ', round(demonstration.gyeonggi.wvs,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(basedata$demonstration[basedata$Region=='Gyeonggi']), main=paste0('Online Sample: Gyeonggi\n(mean = ', round(demonstration.gyeonggi,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))

barplot(table(japan$demonstration[japan$V256 == 'JP: Minami-Kanto region']), main=paste0('WVS: Minami-Kanto Region\n(mean = ', round(demonstration.minamikanto,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(basedata$demonstration[basedata$Region=='Kanagawa']), main=paste0('Online Sample: Kanagawa\n(mean = ', round(demonstration.kanagawa,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(japan$demonstration[japan$V256 == 'JP: Kyushu region']), main=paste0('WVS: Kyushu Region\n(mean = ', round(demonstration.kyushu,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
barplot(table(basedata$demonstration[basedata$Region=='Okinawa']), main=paste0('Online Sample: Okinawa\n(mean = ', round(demonstration.okinawa,1), ')'), ylab='Frequency', names.arg=c('Have\ndone (1)','Might\ndo (2)','Never\ndo (3)'))
dev.off()
